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ABSTRACT 

r 

We  investigate  multi-grid  method*  for  solving  linear  systems  arising  from 
arc-length  continnation  techniques  applied  to  nonlinear  elliptic  eigenvalue 
problems.  We  find  that  the  usual  multi-grid  methods  diverge  in  the 
neighborhood  of  singular  points  of  the  solution  branches.  As  a  result,  the 
continuation  method  is  unable  to  continue  past  a  limit  point  in  the  Bratu 
problem.  This  divergence  is  analysed  and  a  modified  multi-grid  algorithm  has 
been  devised  based  on  this  analysis.  In  principle,  this  nev  multi-grid 

algorithm  converges  for  elliptic  systems  arbitrarily  close  to  singularity  and 
has  been  used  successfully  in  conjunction  with  arc-length  continuation 
procedures  on  the  model  problem.  In  the  worst  situation,  both  the  storage  and 
the  computational  work  are  only  abont  a  factor  of  two  more  than  the  unmodified 
multi-grid  methods.  „ 
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1.  Introduction 

Many  problems  of  computational  interest  can  be  formulated  as 

G(u,X)  =  0,  (1.1) 

where  u  represents  the  'solution*  (i.e.  flow  field,  displacements,  etc.)  and 
X  is  a  real  physical  parameter  (i.e.  Reynold's  number,  load  etc.)  It  is 
required  to  find  the  solution  for  some  X-intervals,  that  is  a  path  of 
solutions,  [u(X),X].  In  this  paper,  we  use  a  class  of  continuation  based  on 
parametrizing  the  solution  branches  by  arc-length,  say  [u(s),X(s)l.  Amain 
advantage  of  these  arc-length  continuation  methods  is  that  most  singular 
points  on  the  solution  branches  can  be  handled  without  much  difficulty. 
Equations  of  the  form  (1.1)  are  called  nonlinear  elliptic  eigenvalue  problems 
if  the  operator  G  with  X  fixed  is  an  elliptic  differential  operator  [21.  For 
nonlinear  elliptic  eigenvalue  problems,  a  major  portion  of  the  computational 
work  in  the  arc-length  continuation  methods  is  spent  in  solving  large  linear 
elliptic  systems.  In  this  paper,  we  investigate  the  use  of  mnlti-grid  [4] 
methods  for  solving  these  linear  systems.  It  turns  out  that  a 

straight-forward  implementation  of  the  multi-grid  methods  fails  in  the 
neighborhood  of  the  singular  points  and  this  usually  prevents  continuation 
past  limit  points.  This  failure  is  analyzed  and  a  modified  multi-grid  method 
based  on  this  analysis  is  devised.  Even  for  very  singular  systems,  the  new 
multi-grid  algorithm  perforates  satisfactorily  and  never  requires  more  than 
about  twice  the  storage  and  computational  work  as  the  unmodified  algorithm. 

The  arc-length  continuation  methods  will  be  described  in  section  2  and 
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the  multi-grid  methods  in  section  3.  In  section  4,  computational  resnltt  for 
a  model  pioblem  are  presented,  together  with  a  description  of  the  difficulties 
encountered  by  the  multi-grid  method  near  a  limit  point.  The  behaviour  of  the 
multi-grid  method  near  singular  points  will  be  analyzed  in  section  5.  The 
modified  multi-grid  algorithms  designed  to  overcome  these  difficulties  are 
described  in  section  6.  The  paper  ends  with  a  summary  in  section  7. 
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2.  Newton* s  Method  end  Continuation  Techniques 

In  this  paper  we  are  concerned  with  aethods  for  computing  a  family  or 

path  of  solutions  of  (1.1).  The  aethods  we  employed  will  be  based  on  some 
version  of  Newton's  method. 

2.1  Newton's  Method 

Given  a  valne  of  X  and  an  initial  guess  u®  for  the  solution  u(X),  we 
perform  the  following  steps  repeatedly  until  1 1 6ux 1 1  <  e  is  satisfied  : 

G*  6uA  «  -  GtuSX)  (2.1) 

u 

ui+1  -  u*  +  Su1.  (2.2) 

In  the  above,  subscripts  denote  partial  derivatives  and  so  G^  denotes  the 
Jacobian  of  the  operator  G  (with  respect  to  u) .  This  procedure  will  generally 
converge  quadratically  when  it  does  converge.  However,  as  is  well  known,  in 
many  instances  it  will  fail  to  converge  when  the  initial  guess  is  not  'close’ 
to  the  true  solution. 

2.2  Natural  Continuation 

A  plausible  procedure  for  overcoming  this  convergence  difficulty  and  also 
for  determining  the  dependence  of  u  on  X  is  to  start  at  a  known  solution 
(u0,Xq)  on  the  solution  curve  and  use  it  as  initisl  guess  for  a  Newton-type 
iteration  to  find  the  solution  for  s  neighboring  point  on  the  solution  curve 
with  X  close  to  Xq.  The  procedure  is  then  repeated.  Ve  can  improve  on  this 
by  computing  the  derivative,  u^,  at  a  known  solution  and  use  it  to  get  a 
better  initial  guess  for  the  nest  value  of  X  in  a  predictor-corrector  fashion. 
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We  call  thi*  a  natural  continuation  procedure  because  it  corresponds  to 
paraaetrizing  the  solution  curve  by  X,  the  naturally  occurring  paraaeter.  A 
specific  fora  of  this  is  the  aore  or  less  well-known: 

Euler-Newton  Continuation  Procedure : 


Given  a  known  solution  (u^.Xq) , 
of  X  as  follows: 

1.  First  coapute  the  derivative  u. 


Gu  nX  =  -  V 


we  coapute  the  solutions  at  nearby  values 


at  (uq.Xq)  froa 


(2.3) 


2.  Perfora  an  Euler  predictor  step: 

u°  =  uQ  +  (X  -  X0>.  (2.4)  \ 

3.  Use  u°  as  initial  guess  in  Newton's  aethod  : 

G1  (ni+1  -  u*)  -  -  G(ui,X)  (2.5) 

u 

until  convergence. 


4.  Use  (u(X),X)  as  the  new  (Qq>Xq)  and  go  back  to  Step  1. 

Note  that  the  coaputation  of  the  derivative  u^  does  not  cause  auch 
coaputational  overhead  because  we  usually  have  the  factorization  of  the 
Jacobian  G&  coaputed  already  in  the  Newton  step.  Using  such  a 
predictor-corrector  aethod  will  often  allow  us  to  take  a  auch  bigger  step  in  X 
and  thus  reduce  the  overall  cost  of  determining  the  dependence  of  u  on  X. 


Unfortunately,  this  procedure  needs  soae  aodification  in  order  to  handle 
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general  nonlinear  systems  became  of  the  possibility  of  existence  of  nonunique 
solutions.  The  nonnniqneness  usually  aianifests  itself  in  the  fora  of 
existence  of  'singular*  points  where  the  Jacobian  Gq  is  singular  (see  Figure 
2-1).  Points  such  as  point  A  in  Figure  2-1  are  called  limit  points  (or 
turning  points)  and  points  such  as  point  6  are  called  bifurcation  points. 
These  singular  points  are  further  characterized  by  the  conditions  that  G.  i 

A 

Range(Gu)  at  a  limit  point  and  that  G^  e  Range (Gu)  at  a  bifurcation 
point  [12] . 

The  difficulties  that  a  natural  continuation  procedure  will  encounter  at 
singular  points  are  three-fold.  First  of  all,  since  Gq  is  singular  at  these 
points,  Newton's  method  will  at  best  be  linearly  convergent,  making  it  much 
more  costly  to  compute  the  solution.  Moreover,  near  a  limit  point,  there  may 
not  exist  a  solution  for  a  given  value  of  X  (tee  Figure  2-2)  and  hence  the 
iterations  must  fail  to  converge.  Lastly,  we  need  some  mechanism  for 
switching  branches  at  a  bifurcation  point. 

2.3  Aro-length  Continuation 

In  the  pseudo  arc-length  continuation  approach  [12],  these  difficulties 
are  overcome  by  not  parametrizing  the  solution  u  by  X.  Instead,  we 
parametrize  the  solution  branches  using  an  arc-length  parameter  s,  and  specify 
how  far  along  the  current  solution  branch  we  want  to  march. 

To  be  more  specific,  we  let  s  be  the  arc-length  parameter,  and  treat  u(s) 
and  X(s)  as  functions  of  s.  We  can  compute  the  'tangent'  [uU^),  X(sq)]  (where 
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the  dots  denote  differentiation  with  respect  to  s)  of  a  known  solution  at  s=sf 


from  the  following  two  equations: 

0, 


Gu  °0  +  X0  Gk 


IIb0M2  +  U0|2  -  1 


0. 


(2.6) 

(2.7) 


Equation  (2.6)  is  obtained  from  differentiating  G(u,X)  =  0  with  respect  to  s 
and  (2.7)  imposes  the  arc-length  condition.  We  could  theoretically  generate 
the  solution  curve  by  integrating  the  initial  value  problem  obtained  by 
solving  (2.6),  (2.7)  for  u(s)  and  X(s).  Although  this  process  is  subject  to 
the  usual  instabilities  inherent  in  solving  initial  value  problems 
approximately,  it  can  be  an  extremely  effective  procedure.  Indeed  our  pseudo 
arc-length  continuation  procedure  can  be  viewed  as  a  method  for  stabilizing 
Euler  integration  of  (2.6),  (2.7). 


In  the  pseudo  arc-length  continuation  procedure,  we  advance  from  sQ  to  s 
along  the  tangent  to  the  solution  branch  and  require  the  new  solution  u(s)  and 
X(s)  to  satisfy 

N(u(s)  ,X(s) )  =  uj(u(s)  -  u(  Sq)  )  +  X0U(s)  -  Us0))  -  (s  -  S{))  =  0.  (2.8) 
In  addition  we  require,  of  course: 


G(u(s) ,X(s))  =  0.  (2.9) 

Equation  (2.8)  is  the  linearization  of  (2.7)  and  as  indicated  forces  the  new 
solution  to  lie  on  a  hyperplane  perpendicular  to  the  tangent  vector  to  the 
solution  curve  at  s^  and  at  a  distance  (s-Sq)  from  it.  Equation  (2.9) 
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Figure  2—3:  Pseudo  Arc-length  Continuation 


solution  curve  on  which 

G(u(o) ,A(o) )  *  0 
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requires  u(s)  end  X(s)  to  lie  on  the  true  solution  curve  (Figure  2-3).  Ve  now 
solve  the  coupled  system  (2.8)  end  (2.9)  for  u(s)  end  X(s),  given  the  step 

size  (s-Sp)  (efficient  stretegies  for  choosing  the  step  size  ere  discussed 
in  [23]).  Ve  use  Newton's  method,  in  which  case  we  have  to  solve  the 
following  linear  system  at  each  iteration: 


Thu  I  Tg 
II  I  U 
A  I  1  =  1 
I  hxl  In! 


6u 

=  - 

6X 

(2.10) 


It  can  be  shown  that  at  limit  points,  where  G  is  singular  and  G.  t 

U  A 

Range(G&),  the  linear  system  in  (2.10)  is  nonsingular  (see  [12])  and  therefore 
Newton's  method  for  the  coupled  system  (2.8)  and  (2.9)  is  well-defined.  Hence 
limit  points  present  no  problem  end  even  quadratic  convergence  is  achievable. 


At  bifurcation  points,  where  Gu  is  singular  end  G^  e  Range(Gu),  things 
are  more  complicated.  In  the  simplest  case  of  only  one  branch  bifurcating 
from  the  main  branch  (simple  bifurcation),  an  additional  higher  order 
condition  involving  G  ,  Gu^  and  G^^  has  to  be  satisfied.  It  can  be 
shown  [12]  that  this  condition,  together  with  (2.6)  and  (2.7)  and  the  left  and 
right  null  vectors  of  Gq,  enable  two  solutions  for  (uq>Aq)  to  be  computed  at  a 
simple  bifurcation  point,  with  one  solution  corresponding  to  each  branch. 
Using  the  appropriate  pair  of  (Up,Xp)  in  (2.8)  allows  branch  switching.  In 
[7]  a  more  detailed  study  of  the  singular  behaviour  and  branch  switching  at 


bifurcation  is  given. 
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In  order  to  solve  the  linear  system  in  (2.10)  by  direct  methods,  several 
approaches  are  possible.  One  wsy  is  to  perform  Gaussian  Elimination  on  the 
inflated  matrix  A,  with  some  form  of  pivoting  to  insure  stability.  But  this 
approach  completely  ignores  the  sparse  structure  which  is  usually  found  in 
G^'s  arising  from  nonlinear  elliptic  eigenvalue  problems.  In  order  to  take 
advantage  of  the  structure  in  G^,  Keller  [12]  suggested  the  following 
block-elimination  procedure: 

Alaorithm  BE:  (block-elimination) 


Solve 

G 

u 

y 

“  G1 

(2.11) 

and 

G 

JX 

z 

=  -  G. 

(2.12) 

Set 

K 

(-N*z-N)/(N,-N*y> 

(2.13) 

and 

6u 

S 

z  -  61  y. 

(2.14) 

Note  that  only  systems  with  the  coefficient  matrix  G^  have  to  be  solved,  so 
structures  in  G^  can  be  exploited.  Moreover,  only  one  factorization  of  Gu  is 
needed.  It  has  been  shown  [27]  that  even  when  Gq  is  becoming  singular. 
Algorithm  BE  produces  iterates  that  converge  quadratically  at  limit  points. 

Continuation  methods  of  various  forms  and  levels  of  sophistication  have 
been  widely  used  in  the  engineering  literature.  For  a  recent  survey  of 
numerical  methods  for  bifurcation  problems,  see  for  example  [18],  The 
approach  taken  here  is  due  to  Keller  [12],  and  has  recently  been  applied  to 
other  problems  in  fluid  mechanics  (  [5],  [6],  [15]  ,  [16],  [25],  [27]  ). 
A  related  approach  suggested  by  Abbott  [1]  corresponds  (in  a  loose  wsy)  to 
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applying  Algoritha  BE  to  the  matrix  A  with  the  laet  column  permuted  into  the 
first  n  columns  so  that  the  corresponding  coefficient  aatrix  in  Equations 
(2.11)  and  (2.12)  becomes  nonsingnlar  even  at  limit  points.  However,  as  has 
already  been  pointed  out,  any  structure  or  symmetry  in  is  lost  in  the 
process,  and  hence  that  approach  seems  unsuitable  for  large  elliptic  systems 
in  two  or  three  dimensions. 


\ 


. . . 
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3.  Multi-Grid  Methods 

3.1  Introduction 

The  class  of  multi-grid  (MG)  nethods  that  we  use  here  is  based  on  work  by 
Bakhvalov  [3],  Brandt  [4].  Federenko  [8].  Hackbush  [10],  and  Nicolaides  [19]. 
We  shall  only  briefly  describe  here  the  particular  MG  algorithms  that  we  have 
used  for  linear  elliptic  problems  that  arise  in  our  treatment  of  nonlinear 
elliptic  eigenvalue  problems. 

The  particular  way  in  which  we  use  the  MG  idea  is  to  use  a  hierarchy  of 

grids,  rather  than  a  single  one,  in  order  to  speed  up  the  convergence  rate  of 

the  solution  process.  The  MG  process  has  some  very  desirable  theoretical 

properties:  for  certain  elliptic  operators  on  an  n  by  n  grid,  it  computes  the 

2 

approximate  solution  to  truncation  error  accuracy  in  0(n  )  arithmetic 
2 

operations  and  0(n  )  storage.  It  seems  natural  to  consider  the  use  of  MG 
methods  for  solving  nonlinear  eigenvalue  problems.  MG  methods  have  been 
applied  to  solution  of  linear  eigenvalue  problems  by  Hackbush  [11]  and 
McCormick  [17]. 

3.2  The  'Cycle  C*  MG  Algorithm 

The  particular  MG  algorithm  that  has  been  used  in  this  study  is  based  on 
the  'Cycle  C'  algorithm  described  in  Brandt  [4].  This  is  an  algorithm  for 
iteratively  solving  the  discrete  equations  approximating  a  linear  elliptic 
problem  on  a  given  grid,  through  interaction  with  a  hierarchy  of  coarser 
grids,  taking  advantage  of  the  fact  that  the  different  discretizations  on  the 
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different  grids  are  all  approximations  to  the  same  continuous  problem.  Ve 
note  that  there  are  other  KG  algorithms  [4]  proposed  for  implementing 
continuation  procedures  outside  of  the  context  of  the  pseudo  arc-length 
framework.  Some  potential  problems  with  these  related  algorithms  are 
discussed  in  section  3.4.  We  do  not  know  how  well  such  MG  algorithms  perform 
and  we  hope  to  carry  out  our  own  investigation  on  such  related  methods  in  the 
future.  In  this  paper.  MG  algorithms  are  uaed  to  solve  the  fine  grid  discrete 
equations  that  arise  in  the  pseudo  arc-length  continuation  procedure. 

Consider  a  hierarchy  of  grids  (G®,Gk,  ....  ,G**),  with  GM  being  the  finest 


one.  defined  on  a  domain  0  with  corresponding  mesh  sizes  (h^  >  hj  >  .  > 

hjj) ,  and  all  approximating  the  same  linear  elliptic  problem  :  | 

LU*  F  on  Q  (3.1) 

U  *  0  on  9C. 

The  discrete  equation  on  a  grid  Gk  is  written  as: 

Lk  Uk  *  Fk  on  Gk  (3.2) 

«=  0  on  30. 


We  are  primarily  interested  in  obtaining  the  approximating  solution  u  on  the 
finest  grid,  and  we  shall  start  with  an  initial  guess  on  G ^  and  apply  a 
standard  relaxation  procedure  such  as  the  Gauss-Seidel  procedure.  It  is  well 
known  that  the  error  is  reduced  rapidly  in  the  first  few  iterations  but  then 
the  reduction  rate  becomes  very  slow.  By  a  frequency  analysis,  it  can  be 
ahown  that  the  fast  reduction  occurs  when  the  residual  (or  the  error)  in  the 


current  iterate  has  large  harmonics  on  the  scale  of  the  grid,  the  so-called 
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high-frequencies.  Now  st  a  stage  in  the  iterative  process  where  the  error 
reduction  rate  slows  down,  let  the  current  iterate  be  uM.  Define  the  error  vM 

in  the  iterate  as  vM  «  D**  -  uM.  Then  the  error  vM  satisfies  the  following 
equation: 


lM  VM  =  f«  _  lM  uM 


rM 

on  G  , 


(3.3) 


v  =  0  on 

The  residual  R**  is  computable  and  hence  the  original  problem  of  solving  for  IT*1 

|| 

can  be  reduced  to  an  equivalent  one  of  solving  (3.3)  for  v  .  There  seems  to 

be  no  obvious  advantage  in  using  (3.3)  rather  than  continuing  with  the 

M  M 

original  relaxation  procedure  with  u  .  However,  if  the  error  v  and  the 
residual  are  smooth  relative  to  (P,  that  is,  if  their  high  frequency 
components  have  been  smoothed  out  by  the  relaxation  procedure,  then  we  can 
approximate  the  solution  of  (3.3)  on  a  coarser  grid,  say  G*1-1,  by  solving  : 


L«-!  v”'1  -  F*'1  -  iJ-V  one”-1, 

K 


v""1  =0 


on 


dGM_1 


(3.4) 


After  this  problem  is  solved  we  can  interpolate  the  solution  v*5-*  onto  GM  to 


get: 


new  uM  *  old  u**  +  wu-i  *m-i  1  *  (3.5) 

where  is  an  interpolation  factor,  normally  taking  the  value  unity  and 

stands  for  some  interpolation  operator  from  G1  to  GJ .  The  solution  process 

U_1 

for  Equation  (3.4)  on  G  usually  costs  considerably  less  than  the  cost  of 
solving  Equation  (3.3)  on  G**.  If  v**  is  indeed  smooth  ^relative  to  G**) ,  then 

11.1  If  ||  |f_l 

G  should  provide  adequate  resolution  for  v  and  hence  v  should  be  a 

|f 

good  approximation  for  v  .  This  principle  of  transferring  to  a  coarser  grid 
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when  convergence  slows  down  can  be  applied  recursively.  Thus  for  example,  we 

U.1 

can  start  with  a  zero  initial  guess  for  v  in  Equation  (3.4)  and  apply  the 
Gauss-Seidel  relaxation  procedure  to  the  iterates  on  G**”1.  When  the 
convergence  slows  down,  we  can  again  transfer  to  the  next  coarser  grid  GU~2, 
and  so  on.  One  can  view  the  whole  process  as  each  grid  smoothing  just  those 
frequencies  in  the  error  that  are  high  relative  to  its  own  mesh  size,  each 
doing  its  job  efficiently  because  these  high  frequencies  are  precisely  those 
that  are  efficiently  smoothed  out  by  relaxation  procedures. 

The  control  of  when  to  transfer  between  grids  can  follow  a  fixed  strategy 
or  an  adaptive  one.  A  fixed  strategy  could  be  of  the  following  kind  (see 
Nicolaides  [19])  perform  p  relaxation  sweeps  on  each  grid  G^  before 
transferring  to  a  coaser  grid  Gk  1,  and  perform  q  relaxation  sweeps  before 
interpolating  back  to  a  finer  grid  Gk+1 .  An  adaptive  strategy  could  be  as 
follow  (see  Brandt  [4])  :  transfer  to  a  coaser  grid  when  the  ratio  of  the 
residual  norm  of  current  iterate  to  the  residual  norm  a  sweep  earlier  is 
greater  than  some  tolerance  q,  and  transfer  to  a  finer  grid  when  the  ratio  of 
the  residual  norm  of  current  iterate  to  the  residual  norm  on  the  next  finer 
grid  is  less  than  another  tolerance  6.  For  simple  problems  like  Poisson's 
equation  on  a  square,  the  overall  MG  efficiency  is  very  insensitive  to  which 
particular  strategy  is  used  and  what  values  are  used  for  (p,q)  or  (q,6).  We 
shall  refer  to  the  above  particular  fixed  strategy  the  (p,q)  strategy  and  the 
adaptive  strategy  the  (q,8)  strategy. 
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3.3  Indefinite  Problems 

In  the  'Cycle  C'  algorithm  jnst  described,  convergence  on  the  lowest 
(coarsest)  grid  G®  is  obtained  by  repeated  relaxation  sweeps.  For  positive 
definite  matrices,  convergence  on  G®  can  be  guaranteed.  For  indefinite 
problems,  however,  convergence  on  G^  cannot  be  obtained  by  repeated  relaxation 
sweeps,  because  the  components  of  the  error  that  correspond  to  eigenfunctions 
with  negative  eigenvalues  will  grow  as  a  result  of  relaxation  sweeps  (see  the 
analysis  in  section  5).  Therefore,  for  indefinite  problems,  a  direct  solution 
(e.g.  Gaussian  Elimination)  must  be  employed  on  the  coarsest  grid.  If  this 
coarsest  grid  is  fine  enough,  it  will  also  provide  corrections  to  those 
growing  components  of  the  iterates  on  all  finer  grids.  However,  too  fine  a 
grid  for  G°  will  increase  the  cost  of  the  direct  solution  procedure.  Hence  a 
little  care  must  be  taken  regarding  the  size  of  the  coarsest  grid  for 
indefinite  problems.  Fortunately,  for  ’not  too  indefinite'  problems,  G°  can 
be  chosen  coarse  enough  so  that  the  direct  solution  on  G®  will  not  affect  the 
overall  efficiency  of  the  NG  procedure  seriously.  Since  indefinite  problems 
occur  frequently  in  nonlinear  elliptic  eigenvalue  problems  and,  in  particular, 
in  our  model  problem,  we  shall  use  such  a  direct  solution  on  G®  whenever 
necessary. 

3.4  Continuation  Methods 

Brandt  [4]  suggested  using  continuation  methods  in  conjunction  with  the 
MG  procedure.  His  main  idea  is  to  use  coarse  grids  for  continuation,  with 
little  work  and  crude  accuracy,  and  only  use  the  finer  grids  at  the  final 


■< 
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continuation  step  to  achieve  higher  accuracy.  We  have  not  pursued  this  idea 
here.  We  believe  that  it  will  work  as  long  as  we  stay  away  from  singular 
points.  Around  a  limit  point,  however,  the  solution  branches  corresponding  to 

different  grids  may  look  like  the  situation  in  Figure  3-1.  If  we  continue  on 

*  * 

the  coarse  grid  to  X  and  try  to  refine  using  the  finer  grid,  while  keeping  X 

fixed,  we  cannot  hope  to  obtain  a  fine  grid  solution  because  X*  is  larger  than 
the  fine  grid  limit  point  X^  (i.e.  no  fine  grid  solution  exists  for  X  >  Xf) . 
In  the  opposite  case,  there  is  no  coarse  grid  solution  at  X  so  we  cannot  get 
started  on  that  grid.  Hence,  in  general,  we  have  to  be  extremely  careful  in 
using  MG  methods  and  continuation  around  singular  points. 
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Figure  3-1:  Linit  Points  for  Different  Grids 
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4.  Application  to  the  Bratn  Problea 

4.1  Brets '<  Problea 

As  e  typical  example  of  an  nonlinear  elliptic  eigenvalue  problem,  we 
consider  the  Bratn  problem  : 

G(u,X)  *  An  +  X  e11  =  0  on  0,  (4.1) 

n  =  0  on  30. 

Equation  (4.1)  arises  in  many  physical  problems,  for  example,  is  chemical 
reactor  theory,  radiative  heat  transfer,  and  in  modelling  the  expansion  of  the 
universe.  The  domain  Q  is  the  nnit  interval  [0,1]  in  R* ,  or  the  unit  square 
[0,l]x[0,l]  in  R^ ,  or  the  unit  cube  [0,l]x[0,l ]x[0,l ]  in  .  There  are  no 
bifurcation  points  in  this  problem,  all  the  singular  points  are  limit  points. 
The  behaviour  of  the  solution  near  the  singular  points  has  been  studied 
numerically  [1,  26]  and  theoretically  [14,  20,  21,  2 43.  Typical  solution 
diagrams  are  shown  in  Figure  4-1.  For  both  the  one  and  two  dimensional  cases, 
the  problem  has  exactly  one  limit  point,  whereas  the  three  dimensional  case 
has  infinitely  many  limit  points  (if  0  is  a  sphere).  From  now  on  we  only 

consider  the  two  dimensional  case,  with  0  the  unit  square.  For  this  case,  the 

value  of  X  and  the  corresponding  !lu|ls  at  the  limit  point  are  given  by  :  X 
=  6.81  and  llull^  =  u(0.5,0.5)  =  1.39.  For  X  >  X*,  Equation  (4.1)  has  no 

solution,  and  for  X  <  X*,  it  has  exactly  two  solutions. 


■  1  -W  .. 


A 
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Figure  4-1:  Solution  for  the  Bretu  Problem 


4.2  Are-length  Continuation  with  Direct  Methods 

We  first  apply  the  arc-length  continuation  method  of  Section  2  to 
Equation  (4.1)  using  direct  methods.  For  this  problem,  a  trivial  solution  is 
(u  *  0,  A  =  0).  We  can  thus  start  at  this  trivial  solution  on  the  lower 
branch  and  march  along  the  solution  branch,  past  the  limit  point,  and  continue 
on  to  the  upper  solution  branch.  Since  the  only  singular  point  in  this 
problem  is  a  limit  point,  this  in  principle  presents  no  problem  to  the 
arc-length  continuation  procedure,  although  the  step  size  might  have  to  be 
reduced  and  controlled  appropriately  near  the  limit  point.  If  desired,  the 
limit  point  can  be  accurately  determined  by  other  related  techniques  [1,  131. 

The  derivatives  of  the  operator  G  in  Equation  (4.1)  that  are  needed  for 
the  arc-length  continuation  technique  are  : 

Gu  -  A  +  k  e*  .  (4.2) 

G^  =  e11  .  (4.3) 

Now  if  we  approximate  the  Laplacian  operator  by  the  standard  five-point 
stencil  on  a  uniform  grid,  the  operator  G^  will  be  approximated  by  the  usual 
block  tridiagonal  matrix  and  the  operator  G^  by  a  column  vector. 

In  the  application  of  the  arc-length  continuation  technique,  we  will  have 
to  repeatedly  solve  linear  systems  of  equations  with  the  matrix  given  by  G^. 
The  solution  of  these  linear  systems  is  the  central  part  of  the  arc-length 
continuation  method.  Bence,  an  efficient  linear  system  solver  is  crucial  to 
the  overall  performance  of  the  continuation  technique.  In  this  section,  we 
present  some  computational  results  for  Bratu's  Problem  using  a  direct  method 
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(Gaussian  Elimination)  of  solution  of  the  linearized  difference  equations. 
For  large  problems,  this  would  be  prohibitively  expensive.  However,  the 
results  here  are  intended  to  demonstrate  the  performance  of  the  continuation 
procedure  independent  of  the  linear  algebra  method  employed.  In  the  next 
section,  we  shall  investigate  the  use  of  multi-grid  methods  for  solving  the 
linear  equations.  It  should  be  pointed  out  that  is  generally  not 
separable,  and  therefore  we  cannot  use  fast  Poisson  solvers  directly  even  on 
rectangular  domains.  Moreover,  this  matrix  is  indefinite  on  the  upper  branch, 
and  hence  iterative  methods  like  Successive-Over— Relaxation  cannot  be  used 
directly. 

We  present  some  of  our  computed  results  in  Figure  4-2.  Only  the 
behaviour  of  the  solution  branch  near  the  limit  point  for  a  few  relatively 
coarse  discretizations  are  presented.  These  are  to  be  compared  with  the 
values  :  X*  =  6.80811698  and  u(.5,.5)  =  1.3916603  for  a  grid  with  h  *  1/24 
with  the  nine-point  finite  difference  operator  as  computed  by  Abbott  [1]  and 
to  the  easily  obtainable  exact  solution  (X  *  18/e  *  6.62183,  u  =1)  for  the 
case  h  =  1/3.  As  expected,  the  step  size  3s  *  s  -  sQ  had  to  be  suitably 
controlled  near  the  limit  point,  but  otherwise  we  encountered  no  difficulty  in 
continuing  past  the  limit  point. 

4.3  Arc-length  Continuation  with  Multi-Grid  Methods 

In  this  section  we  discuss  the  use  of  MG  methods,  rather  than  direct 
methods,  for  solving  the  linear  equations  that  arise  in  the  continuation 
procedure.  The  MG  method  that  we  use  wss  described  in  Section  3  and 


Figure  4-2:  Computed  Results  for  Brstu’s  Problem  Near  Limit  Point 
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Gauss-Seidel  is  the  smoothing  relaxation  process.  Since  the  Jacobian  aiatrix 
Gu  becomes  indefinite  on  the  spper  branch,  we  ose  a  direct  method  on  the 
coarsest  grid  in  the  neighborhood  of  the  the  limit  point  and  on  the  upper 
branch. 

We  started  the  continnation  procedure  with  the  trivial  solution  (u  “O,  X 
=  0),  with  h  =  1/4  on  the  coarsest  grid,  and  a  total  of  four  levels  of  grids, 
making  the  finest  grid  with  h  *  1/32.  As  expected,  the  MG  method  worked  fine 
and  we  were  able  to  continue  up  to  very  close  to  the  limit  point,  at  X  *  6.804 
on  the  lower  branch.  However,  we  noticed  that  the  convergence  of  the  MG 
method  deteriorates  as  we  move  in  towards  the  limit  point.  For  example,  the 
number  of  equivalent  relaxation  sweeps  on  the  finest  grid  required  to  reduce 
the  residual  norm  by  an  order  of  magnitude,  which  is  a  convenient  way  of 
measuring  the  efficiency  of  MG  methods,  went  from  about  5  at  X  *  0  to  about  20 
at  X  *  6.803  and  to  divergence  at  X  -  6.805.  The  divergence  occurred  in  the 
MG  method  and  not  in  the  Newton  iteration.  It  is  not  due  to  the  possible 
indefiniteness  of  the  Jacobian  matrix  on  the  finest  grid.  This  can  occur  near 
the  limit  point  after  a  large  Euler-predictor  step.  But  we  performed  other 
tests  starting  on  the  upper  branch,  away  from  the  limit  point,  where  the 
Jacobian  matrix  is  indefinite,  and  the  MG  method  performed  as  efficiently  as 
on  the  lower  branch.  From  our  experience,  this  divergence  is  strictly  a 
phenomenon  associated  with  the  limit  point,  and  to  the  best  of  our  knowledge, 
has  never  been  discussed  or  analysed  in  the  literature,  le  study  this  effect 


in  section  5 . 


The  exact  valae  of  X  at  vhich  thia  divergence  first  occurs  varies 
slightly  with  the  size  of  the  coarsest  grid  hQ,  but  is  quite  independent  of 
the  other  paraaieters  of  the  Cycle  C  algorithm  (e.g.  tj  and  6).  In  all  the 
cases  we  have  run,  this  divergence  made  it  iapossible  to  continue  past  the 
limit  point.  Therefore,  a  remedy  is  needed.  Before  we  can  find  one,  we  must 
understand  the  reason  for  the  divergence. 


! 
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5.  Analysis  of  Multi-Grid  Methods  for  Near-singular  Systems 

For  the  present  analysis,  we  as sane  that  the  linear  operator  L  is 
self-adjoint  and  has  the  complete  set  of  orthonroaal  eigenfunctions  (g, 

....  }  with  corresponding  real  eigenvalues  (pj  <.  ^2 The  operator  Gu  in 
the  Bratn  problem  clearly  satisfies  the  above  hypothesis.  Thus  the  solution  0 
to  L  D  *  F  can  be  written  as: 

CD 

XJ  ■  ^  a.  aj  -  <^.F>,  j-1,2,...  (5.1) 

We  assnae  that  the  discrete  approximations  Lk  to  the  continuous  L  are 
symmetric.  Thus  they  have  real  eigenvalues  (pk  <  pk  <.  ...  <  pjj)  and  a 
complete  set  of  orthonormal  eigenvectors  {£k,  $k,  •••  •  •  Here  is  the 

dimension  of  the  matrix  representing  Lk.  For  most  reasonable  approximations, 
and  certainly  for  the  five  point  formula  used  for  the  Bratu  Problem  on  a 
rectangle  this  is  true. 


Assume  that  after  iterating  (relaxing)  on  the  grid  Gk,  convergence  has 
slowed  down  and  a  transfer  to  the  next  coarser  grid  is  desired.  Let  the 
current  iterate  be  uk,  and  the  corresponding  'correction'  be  vk  so  that  U*  * 
uk  +  vk  where  Uk  satisfies  LkUk  «  F*.  The  correction  problem  is  given  (as  in 
section  3)  by: 

Lk  vk  ■  Kk  =  Fk  -  Lk  uk,  in  Gk;  vk  *  0  on  8Gk.  (5.2) 

This  is  approximated  on  Gk  1  by 

Lk-1  vk-l  .  jk^  Rk  in  Gk.  vk-l  =  o  on  9Gk_1.  (5.3) 

Using  the  eigenvector  expansion  of  vk  in  (5.2)  we  get: 


•  t  ■' 


(5.4) 


-  28  - 


where 

ai  "  <**'**>  1  (5.5) 
Suppose  now  that  (5.3)  is  solved  exactly  (by  either  direct  solution  or  Cycle  C 
or  any  other  means)  on  Gk~*.  The  solution  vk~*  is  then: 


where 


N, 


„k-l 


1 


k-1 


k-1 


(5.6) 


1  -  <lJ"1Rk,?^_1>  /  Uk_1. 


The  key  idea  in  the  MG  method  is  that  if  vk  and  Rk 


,k-l 


Thus  it  is 


be  well  approxiaiated  on  G‘ 

considerations  that 

Tk  k-1  -  k 

Vl  v 

Using  (5.4)  and  (5.6),  this  is  equivalent  to: 
Nr. 


k-1  Tk 


xi-i  h 


k-1  ~ 


k  *k 


li  -i 


This  will  be  the  case  if 

.k-1  -  «k 


U)  ri-i  «i  “  ««Vl  • 

UiiVl' 


(b) 


(5.7) 

are  smooth  enough,  they  can 
important  for  efficiency 

(5.8) 

(5.9) 

(5.10) 

(5.11) 


3  ~ 

We  shall  use  the  symbol  to  mean  rather  loosely  'approximately  eqnal 
to’.  The  meaning'' should  be  clear  by  context.  Also,  we  shall  assume  that  the 
interpolation  fafctor  w.  .  in  Equation  (3.5)  is  equal  to  one  unless  stated 
otherwise. 
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(c)  ak  =  0.  i  >  Nk_j.  (5.12) 
Conditions  (5.10)  and  (5.11)  ensure  that  the  coarse  grid  correction  vk_* 
improves  the  lower  modes  of  the  iterate  uk.  Condition  (5.12)  is  essentially 
the  smoothness  required  of  vk  on  Gk  (i.e.  negligible  higher  modes). 


Now  condition  (5.10)  is  satisfied  for  the  low  frequency  eigenfunctions  of 

k  v—  i 

the  continuous  operator  L  if  the  grids  G  and  G  are  both  fine  enough  to 

resolve  these  eigenfunctions.  This  holds  in  many  cases  since  the  lower 

eigenfunctions  of  most  second  order  elliptic  operators  over  smooth  domains  are 
very  smooth.  For  the  Bratu  problem,  the  eigenfunctions  are  very  close  to 
products  of  sines  and  cosines  (the  eigenfunctions  of  the  Laplacian  operator) 

and  so  the  lower  modes  are  easily  resolved  by  very  coarse  grids.  Condition 

(5.11),  on  the  other  hand,  turns  out  to  be  violated  if  the  operator  L  is  near 
singular.  This  is  what  caused  the  divergence  of  the  Cycle  C  algorithm  in  the 
arc-length  continuation  procedure  as  we  approach  the  limit  point  (see  section 
4.3).  We  shall  analyse  this  case  next. 


From  (5.5)  and  (5.7),  condition  (5.11)  becomes: 
<Ik_1Rk,{k_1>  /  pk_1  -  <Rk,?k>  /  pk, 
i  1  i  i  Nk-1 . 


(5.13) 


We  dais  that  if  condition  (5.10)  is  satisfied,  and  if  the  transfer  fron  G*  to 
G^~*  is  done  only  after  the  residual  Rk  has  been  smoothed,  then  the  numerators 


in  (5.13)  will  have  approximately  the  saae  value*  To  show  this,  we  expand  R1 


as 
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(5.14) 


where 

r£  -  <Rk,?k>.  (5.15) 
Thus  the  numerator  on  the  right  hand  side  of  (5.13)  is  precisely  r ^ ,  To 
estiaiate  the  numerator  on  the  left  hand  side  of  (5.13).  we  proceed  as  follows: 


rk-l  rk 


rk-1 


i  ‘k 


'Nk-1+1 


(5.16) 


Now  if  condition  (5.10)  holds,  its  converse: 

Ik"1  ?i  "  ?i_1 '  1-^i-^Nk-l '  (5-17) 

V  V 

also  holds.  Also,  if  R  has  been  smoothed  on  G  ,  then  ri  [for  Nk_j<i.£Nk]  must 
be  small  compared  with  tj  [for  l.£i<Nk_j]  ).  Alternatively  (5.12)  assumes  ak  ■ 
r^p.  *  °  *or  i  >  ^k-l  ‘  Tko^ofore,  we  can  approximate  in  (5.16)  by  dropping 
the  second  sum  on  the  right  hand  side  to  get 


Hence 


(5.18) 


<Ik_1  Rk,  ?k_1>  £  r.,  lii^.j.  (5.19) 

Therefore,  from  (5.15)  and  (5.19),  we  have,  as  claimed  earlier, 

<lk_1Rk.?k~1>  £  <Rk.?k>  for  liil^.j.  (5.20) 


The  relations  in  (5.20)  imply  that  condition  (5.13)  will  be  true  if 
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M*  /  Mi""1  =  1.  liiiVj.  (5.21) 

Actually,  these  conditions  need  to  be  strengthened  in  order  to  guarantee  that 

V_1  V 

the  visit  to  G  actually  improves  the  accuracy  of  u  .  This  can  be  seen  as 

v  k— t 

follows.  The  error  in  the  iterate  u  before  the  transfer  to  G  is  given  by 


"k 

k  \  k  k 
old  error  ■  v  *  >  *T 

H 


(5.22) 


K  t— 1 

From  (3.5),  the  new  error  in  u  after  coming  back  from  a  visit  to  G  i 
given  by 


is 


rk  _k-l 
"k-1  *k-l 


(5.23) 


In  view  of  (5.4)  and  (5.6),  the  above  gives: 
N». 

new  error  ~ 


=  2  (a^-w^jS*  J)?i  +  higher  modes 

=2  <1  -  «i-1  /  *i>  *i  ?i  +  h igher  modes 

From  (5.5),  (5.7)  and  (5.20),  we  have 

•r  /  *s :  *5  /  -r- 

and  therefore  we  can  write  the  new  error  in  (5.24)  as: 


(5.24) 


new  error 


'  ^  (1  "  wk-l>1i/»‘i'1),i5i  + 

i*l 


(5.25) 


higher  modes. 

For  obvious  efficiency  and  convergence  considerations,  the  new  error  should 
preferably  be  less  than  the  old  error,  at  least  for  the  lower  modes.  In  other 


words,  condition  (5.21)  should  be  strengthened  to 

k,..k-li 


i.e. 


11  ■  ViY"!  1  <  l* 


o  <  Vi  ^i^i"1  <  2* 


for  l<ilNk_1 , 


(5.26) 


(5.27) 
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Now  if  the  ratio*  of  eigenvalues  in  (5.21)  are  not  close  to  unity,  the 
interpolation  factors,  w^_j#  should  be  chosen  so  that  condition  (5.27)  is 
satisfied.  Otherwise  the  new  error  can  be  larger  than  the  old  error  in  tone 
modes. 

It  should  be  pointed  out  that,  in  general,  condition  (5.27)  is  not 
necessary  for  the  convergence  of  the  Cycle  C  algorithm.  This  is  the  case,  for 
instance,  if  L  and  the  L*'s  are  all  positive  definite.  Then  Gauss-Seidel 
sweeps  on  any  grid  Gk  will  reduce  the  amplitude  of  every  mode  present  in  the 
error.  In  such  cases,  convergence  on  any  grid  can  be  achieved  by  merely  doing 
enough  relaxation  sweeps.  Then  it  is  not  necessary  for  the  next  coarser  grid 
to  provide  any  improvement  on  the  current  iterate,  although  it  would  obviously 
improve  the  efficiency  of  the  overall  algorithm  if  it  does  so.  In  fact,  the 
MG  method  derives  its  efficiency  from  the  very  fact  that  the  coarser  grids  do 
provide  improvements  in  the  current  iterate  uk  in  the  lower  modes.  These  are 
precisely  those  modes  that  have  poor  convergence  rates  for  the  relaxation 
sweeps  on  Gk.  Thus,  even  in  the  positive  definite  case,  it  is  important  (from 
an  efficiency  viewpoint)  that  conditions  (5.27)  hold,  at  least  for  small  i’s. 

If  the  operator  L  and  the  Lk's  are  indefinite  the  situation  is  different 
because  some  modes  will  grow  if  we  simply  perform  relaxation  sweeps  on  a  fixed 
grid.  Such  modes  have  to  be  corrected  by  going  to  coarser  grids  and  using  a 
direct  method  on  the  coarsest  grid.  Further  the  interpolation  factors,  w^_^, 
should  be  chosen  such  that  condition  (5.27)  is  satisfied  for  these  modes. 
Condition  (5.27)  has  been  suggested  by  Brandt  [4]  for  indefinite  problems. 
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However  as  we  show  later,  most  nonlinear  eigenvalue  problems  with  limit  points 
and  bifurcation  points  abound  with  indefinite  operators  but  they  do  not  cause 
difficulties  in  the  sense  of  violating  condition  (5.27).  Essentially  only  one 
mode  causes  problems  on  each  Gb  and  it  is  the  mode  that  correspond  to  the 
eigenvalue  that  is  nearest  zero  as  the  singular  point  is  approached.  Merely 
including  the  interpolation  factors  so  that  condition  (5.27)  is  satisfied 
turns  out  to  be  very  inefficient.  Further,  it  is  not  clear  that  such  factors, 
wk-l'  c,n  be  f°und  at  all  in  this  case. 


Another  source  of  difficulty  is  that  the  process  of  interpolating  vk_1 
into  Gb  introduces  high  frequency  errors.  That  is,  the  exact  relation 
corresponding  to  (5.10)  is: 


tk  k-i  =  Ek  r  k  .k 

rk-i  -  «i +  2.  bij  * 

J 


i  ■=  1.2, 


,N. 


k-1  ’ 


for  lliiN^j,  (5.28) 

V 

and  the  coefficients  b„  may  be  large  for  j  >  N^_j .  This  would  result  in  a 
violation  of  (5.12).  Fortunately,  these  high  frequency  errors  are  very 
efficiently  smoothed  out  by  the  subsequent  relaxation  sweeps  on  Gb,  and  thus 
these  errors  are  automatically  corrected. 


For  elliptic  operators  which  are  ’far’  from  being  singular  and  with  a 
reasonable  grid  system  (GkJ  condition  (5.27)  can  be  assured.  For  example,  if 
L  is  the  negative  Laplacian,  -A,  on  a  unit  square  with  Dirichlet  boundary 
conditions,  then  it  is  known  (e.g.  [9])  that  the  eigenvalues  of  L  are  given 
by 


1 
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M_  _  =  (but)2  +  (an)2.  (5.29) 

id  » n 

The  corresponding  eigenfunctions  are: 

^m.n  =  sin  sin  •  (5.30) 

These  eigenfunctions  evaluated  at  the  discrete  interior  grid  points  of  a 
uniform  mesh  on  the  unit  square,  give  the  eigenfunctions  of  the  discrete 
5— point  approximations,  ■  _^h*  with  h  being  the  uniform  mesh  size.  The 

eigenvalues  of  Lk  are,  with  6x  =  6y  =  h^: 

Mm  n  =  4fsin2(®nllk/2)  +  sin2<nnhk/2)J  /  h£.  (5.31) 

Some  of  these  eigenvalues  are  tabulated  in  Table  5-1  for  various  mesh  sizes, 
k  k— 1 

h.  .  The  ratios  p  / u  are  given  in  Table  5-2.  We  see  from  Table  5-2  that 
£  Di&  n,ii 

condition  (5.27)  is  satisfied,  with  w^  =  1,  for  all  lower  modes  shown. 

These  ratios  are  very  close  to  unity,  even  for  the  case  where  the  coarsest 

grid  has  only  one  interior  point.  We  have  seen  from  condition  (5.11)  that 

this  closeness  to  unity  is  very  desirable  and  this  fact  partly  explains  the 
well-documented  success  of  MG  methods  for  the  Laplacian  operator. 

Near  the  limit  point  of  the  Bratu  problem,  the  operator  L  *  =  a  +  keB 

behaves  very  much  like  a  shifted  Laplacian  operator.  Clearly,  if  the  factor 

eB  were  replaced  by  a  constant,  a  say,  then  G^  is  replaced  by  the  the 
Laplacian  operator  with  a  shift  aX.  Then  the  eigenvalue  ratio  p^ 
valid  for  aX  =  0,  is  replaced  by: 

(41,1  "  “  «*>•  (5.32) 

Since  0  <  u  <  1.4  the  factor  en  does  not  vary  much  and  we  assume  this 


approximation  to  be  valid  for  some  a  >  0. 


The  situation  is  depicted 
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Table  5-1: 
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Table  5-2: 


Ratios 


k-1 

0,0 
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Figure  5-1:  Spectrum  of  Shifted  Lepleciea 
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graphically  in  Figure  5-1  for  the  grid  system  that  was  used  for  Table  5-1.  As 
the  shift  aX  approaches  the  group  of  eigenvalues  corresponding  to  the  (1.1) 
siode  from  below,  the  ratios  in  (5.31)  increase.  As  aX  continues  to  increase 
the  ratio  of  eigenvalues  will  become  greater  than  2,  then  increase  towards  +«, 
jump  to  -®  discontinuosly ,  and  start  increasing  from  -®  to  1.  The  situation 
is  depicted  in  Figure  5-2. 

Ve  thus  see.  under  the  above  assumptions,  that  condition  (5.27)  is  first 
violated  by  the  lowest  mode  (i.e.  the  (1,1)  mode)  on  the  two  coarsest  grids  G® 
and  G* .  In  fact  the  lowest  eigenvalues  for  the  Bratn  problem  computed  at  the 
first  point  on  the  solution  branch  where  Cycle  C  diverged,  yields  the  ratio 
almost  exactly  2!  On  the  other  hand,  even  at  this  point,  condition  (5.27)  is 
satisfied  by  the  (1,1)  modes  on  the  finer  grids.  In  other  words,  the 
divergence  of  Cycle  C  is  seen  to  be  caused  by  one  near-singular  grid  out  of 
the  whole  hierachy  of  grids  present.  The  mode  that  becomes  singular  at  the 
limit  point  of  the  Bratu  problem  is  the  (1,1)  mode,  and  this  occurs  first  on 
the  G®  grid.  As  the  limit  point  is  approached,  Lk  on  some  of  these  grids  may 
even  become  indefinite,  while  others  (the  finer  grids)  may  still  be  positive 
definite.  Essentially,  the  near-singular  grid  causes  the  (1,1)  mode  component 
of  the  correction  vk  *,  when  viewed  as  an  approximation  to  vk,  to  have  the 
right  direction,  but  the  wrong  magnitude .  This  phenomenon  is  not  limitted  to 
the  Bratu  problem.  The  only  thing  special  about  this  problem  is  that  it  is 
the  eigenvalue  of  the  (1,1)  mode  that  becomes  zero  at  the  limit  point.  For 
other  problems,  the  eigenvalue  of  the  operator  L  that  becomes  zero  as  the 
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singular  point  is  approached  night  correspond  to  other  nodes.  Although  the 
singular  point  in  the  Bratn  problem  is  a  linit  point,  ve  can  expect  the  sane 
behaviour  at  a  bifurcation  point. 

Having  now  understood  the  cause  of  the  divergence  of  the  HG  nethod,  in 
the  next  section  we  shall  discuss  some  modifications  to  the  basic  Cycle  C 
algorithm  that  are  designed  to  overcone  such  difficulties. 
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6.  Kntdiu  sad  New  Algor ithas 

la  this  section  we  discuss  approaches  that  have  been  devised  to  overcoae 
the  difficulties  with  the  MG  method  near  singular  points.  The  first  goal  is 
to  modify  the  basic  Cycle  C  algorithm  so  that  it  will  converge  for  values  of  X 
close  enough  to  the  liait  point  so  that  the  arc-length  continuation  procedure 
can  take  us  past  the  liait  point  onto  the  upper  solution  branch.  A  more 
ambitious  goal  is  to  aodify  Cycle  C  further  so  that  it  will  converge 
arbitrarily  close  to  the  singular  point.  Such  an  algorithm,  when  used  in 
conjunction  with  the  arc-length  continuation  technique  for  tracing  solution 
branches,  will  Bake  the  overall  algorithm  much  more  robust.  Moreover,  such  an 
algorithm  may  prove  to  be  useful  for  locating  singular  points  accurately, 
either  using  an  arc-length  continuation  based  procedure  [13],  or  some  other 
procedure  that  uses  the  operator  near  the  singular  point  [22],  We  shall 
see  that  the  first  goal  is  relatively  easy  to  achieve,  whereas  the  second  goal 
is  much  more  difficult.  However,  we  have  devised  a  Cycle  C  based  algorithm 
that  has  performed  very  well  when  applied  very  close  to  the  limit  point.  The 
approaches  that  we  have  tried  and  that  lead  to  the  final  algorithm  will  be 
discussed  in  this  section.  We  shall  describe  them  in  the  sequence  that  they 
were  tried. 

Before  we  proceed,  however,  we  have  to  explain  a  few  general  strategies 
that  were  used.  First  of  all,  Gauss-Seidel  and  many  other  relaxation  schemes 
are  not  very  effective  in  smoothing  the  lower  modes,  especially  modes  with 
near  zero  eigenvalues.  Hence,  these  modes  must  be  eliminated  by  means  other 


than  relaxation,  even  on  the  coarsest  grid.  Therefore,  unless  stated 
otherwise,  we  shall  use  a  direct  solution  on  the  coarsest  grid  even  though  the 

operators  L  's  nay  be  positive  definite.  This  does  not  affect  the  overall 
efficiency  very  much  because  the  coarsest  grid  has  so  few  points  that  direct 
solution  is  very  fast  and  efficient. 

Another  strategy  concerns  the  treatment  of  the  mode  that  causes  the 
divergence;  that  is  the  mode  with  a  near  zero  eigenvalue,  say  ^ .  in  all  the 
algorithms  that  are  discussed,  this  mode  is  treated  separately  from  the  other 
modes.  To  do  this,  it  is  essential  to  have  approximations  to  this  mode  and  to 
its  corresponding  eigenvalues,  say  5^  and  |ij,  respectively.  Here  we  have  to 
strike  a  balance  between  accuracy  and  efficiency.  If  we  compute  the 
exactly,  then  we  can  completely  eliminate  the  error  components  on  each 
grid.  Thus,  the  problem  on  G  can  be  reduced  to  one  in  which  is  zero  (see 
(5.25)).  When  this  is  done,  we  do  not  need  to  satisfy  condition  (5.27)  for 
this  mode.  On  the  other  hand,  the  work  involved  in  computing  accurate 
approximations  to  p*  and  for  each  k  would  be  at  least  as  much  as  solving 
the  original  linear  system.  Our  compromise  has  been  to  compute  an 
approximation  5^  to  on  the  coarsest  grid,  G®,  by  a  few  steps  of  inverse 

iteration  with  zero  shift  (since  the  eigenvalue  we  want  is  near  zero).  This 
is  very  inexpensive  since  G®  is  quite  coarse  and  the  LU  factors  of  L®  are 
already  available.  Then  we  interpolate  onto  the  finer  grids.  To  eliminate 
the  high  frequency  errors  introduced  in  these  interpolations,  we  do  two 
things:  1)  use  higher  order  interpolation,  e.g.  cubic  instead  of  linear,  2) 
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smooth  the  interpolated  eigenfunctions  by  performing  a  few  relaxation  sweeps 
on  {j  «  0.  Estimates  of  the  eigenvalues,  p*.  are  then  computed  using  the 
Rayleigh  Quotients:  <$*,L*£*>.  Ve  view  this  as  a  preprocessing  phase  of  the 
algorithm  and  the  extra  work  is  usually  minimal  compared  to  the  overall  work. 
Furthermore,  since  the  eigenfunctions  (not  the  eigenvalues)  do  not  change  very 
much  in  the  neighborhood  of  the  singular  points,  we  can  use  the  same 
approximation  for  different  linearized  operators  L*.  The  storage  required  to 
store  these  eigenfunctions  is  less  than  twice  the  size  of  the  finest  grid. 

Ve  use  the  (q,5)  adaptive  version  of  the  Cycle  C  algorithm,  unless 
otherwise  stated.  The  first  modified  algorithm  is  the  following. 

6.1  Under-  and  Over-  Interpolation 

The  idea  is  to  choose  wk-J  in  (3.5)  for  interpolation  onto  Gk,  such  that 
condition  (5.27)  is  satisfied  for  ,  Clearly  the  value: 

wk-1  *  (Ij  1/Jij  .  (6.1) 

is  in  some  sense  optimal  since  it  eliminates  the  term  in  (5.25).  For  the 
case  discussed  in  Section  4.3,  this  modification  allows  the  computation  to 
continue  past  the  point  X.  ■  6.804,  where  divergence  of  Cycle  C  first  occurred. 
In  fact  (with  a  little  luck)  we  succeeded  in  continuing  around  the  limit  point 
onto  the  upper  branch.  Sere  the  eigenfunction  no  longer  presented 

difficulties  for  the  MG  algorithm.  For  some  of  these  cases  p®  is  actually 
negative  and  therefore  (6.1)  yields  a  a  negative  value  for  w1 .  In  this  case 
the  transfer  from  G°  to  G1  violates  condition  (5.27)  for  all  modes  other  than 
(j.  The  errors  in  these  modes  must  be  reduced  by  extra  relaxation  sweeps  on 


g 


G1 .  In  other  words  G®  only  provides  s  proper  correction  on  G*  for  the 
■ode,  >11  higher  nodes  are  treated  incorrectly  during  the  transfer.  The 
efficiency  of  the  algorithn  thns  suffers.  This  effect  is  especially 
pronounced  if  sone  factors  w^  are  either  very  large  or  negative  or  (worse) 

both.  The  algorithn  is  very  sensitive  to  the  paraneters  ( rj , 6 )  and  thns  is  not 

robust.  It  can  even  diverge  if  the  higher  nodes  are  not  reduced  fast  enough 
on  G*  after  the  transfer  fron  G*~*. 

Even  worse,  the  above  algorithn  will  not  work  for  indefinite  problens  in 
which  sone  intemediate  eigenvalue  is  near  zero.  For  exanple,  if  the  spectra 
of  the  Lk  are  sinilar  to  those  in  Figure  6-1,  the  interpolation  factors  w^  are 

controlled  by  the  £k  belonging  to  eigenvalues  fik  near  zero.  On  the  other 

hand,  the  eigenfunctions  {^require  that  condition  (5.27)  be  satisfied  because 
these  nodes  cannot  be  liquidated  by  relaxation.  Conflicts  can  occur  when  £k 
requires  to  be  negative  while  {^requires  w^  to  be  positive.  Indefinite 
problens  of  this  type  occur  frequently  in  nonlinear  eigenvalue  problens.  Mere 
under-  or  over-interpolation  nust  run  into  difficulties  for  such  problens, 
near  the  singular  points. 

The  above  considerations  nake  it  clear  that  the  eigenfunction  with  the 
near-zero  eigenvalue  nust  be  isolated  and  treated  different  fron  the  other 
eigenfunctions.  We  use  the  approxinate  eigenfunctions  that  are  conputed  in 
the  preprocessing  phase  for  this  purpose  in  the  following  procedure. 
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6.2  Under-  and  Over-  Interpolate  the  Singular  Eigenfunction  Only 

Ve  use  an  interpolation  different  from  that  in  (3.5).  Specifically  if 


k-1  V 

on  6  ,  we  interpolate  it  onto  G  by 


k  _  „  k-lTk  ,k-l  .  Tk  V  k-l.k-1 
*k-l *1  Ik-l**l  Xk-l2  „  #i  5i 


(6.2) 


(6.3) 


Further  i*  choaen  to  satiafy  (6.1).  Since  we  only  have  an  approximation 

to  £*,  we  use,  inatead  of  (6.3): 


vk  "  I**  tv*'1  - 


(6.4) 


In  practice,  this  performed  much  better  than  indiscriminate  under-  and 
over- interpolation  described  in  section  6.1.  It  was  the  more  efficient  when 
both  procedures  worked.  In  many  cases  when  (6.1)  yields  large  and/or  negative 
values  for  w^,  only  the  current  scheme  converges.  In  principle,  it  will  also 
work  for  indefinite  problems  like  that  depicted  in  Figure  6-1.  The  efficiency 
in  most  cases  was  very  respectable;  in  the  range  of  6-10  units  per  order  of 
magnitude  reduction  in  the  residual.  It  is  also  quite  insensitive  to  the 
parameters  (q,6).  Thus,  it  can  be  nsed  very  efficiently  and  reliably  with  the 
arc-length  continuation  procedure  for  tracing  out  solution  branches. 

Unfortunately,  this  improved  algorithm  fails  when  the  magnitude  of  w^ 
becomes  too  large.  This  occurs  when  L*  is  very  nearly  singular,  that  is  with 
pk  very  close  to  zero.  Since  we  only  have  an  approximation  to  $*,  large 
factors  w^  in  (6.4)  introduce  very  large  errors  in  the  other  modes.  Moreover, 
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the  estimates  |i*  using  Rayleigh-Quotients  tend  to  be  too  large  (relatively) 
when  p*  is  very  snail.  Then  (6.1)  gives  a  value  of  w^  that  is  too  snail. 
Both  of  the  above  result  in  lower  efficiency  and  reliability.  In  extreme 
cases,  this  makes  the  algorithm  impractical.  To  overcome  this  difficulty,  we 
devise  an  algorithm  that  will  work  even  if  one  of  the  operators  L*  is  very 
nearly  singular.  For  this  we  employ  the  idea  of  skipping  a  grid. 

6.3  Skipping  the  Singular  Grid 

The  previous  algorithm  fails  if  the  operator  is  very  nearly  singular  on 
one  of  the  grids,  say  G^.  The  idea  here  is  to  simply  delete  this  grid  from 
the  hierachy  of  grids  used  by  the  MG  algorithm.  If  the  remaining  grids  are 
not  as  singular  as  the  deleted  grid  it  would  seem  that  the  algorithm  described 
in  6.2  should  work.  However,  calculations  show  that  skipping  a  grid  can  cause 
other  problems.  When  G  is  skipped,  the  mesh  changes  more  drastically  from 
Gk  1  to  G*+1,  and  hence  the  interpolation  in  (6.4)  (now  I***  instead  of  l£_j) 
introduces  larger  errors  into  the  higher  modes  on  Gk+1 .  These  high  frequency 
errors  can  cause  divergence  of  the  MG  process  unless  controlled  properly  by 
the  parameters  (q,6) .  A  large  value  of  q,  say  between  .8  and  .9,  makes  the 
algorithm  more  robust  but  involves  more  work  than  for  a  smaller  value  of  q, 
say  .5.  We  encountered  a  case  where,  with  all  else  the  same,  the  new  skipping 
algorithm  converges  for  q  c  .9  but  diverges  for  q  *  .6.  Granted  with  q  *  .9 
the  algorithm  may  be  very  reliable,  such  sensitivity  to  one  parameter  is  very 
undesirable.  Therefore,  we  considered  the  following  modification. 


-  48  - 


6.4  Skipping  the  Singular  Grid  for  the  Singular  Eigenfunction  Only 

The  idea  is  to  skip  the  singular  grid  Gk  for  ^  only,  and  to  keep  it  for 
smoothing  the  other  nodes.  In  the  actual  inplenentation,  we  modify  the 
algorithm  described  in  section  6.2  to  use 

*k_!  «  i^'1  /  4^+1  (6.5) 

f°r  and  *  1  for  all  other  modes  to  transfer  from  Gk  1  to  Gk  and,  after 

k  k+1 

a  few  smoothing  sweeps  on  G  ,  transfer  to  G  with  »k  =  1  for  all  modes. 

k  k 

Note  that  we  do  not  try  to  solve  the  6  equations  for  v  .  Trying  to  do  that 

k  k  k 

would  result  in  large  magnification  of  the  4*  component  in  v  ,  since  is 

k+1 

near  zero.  This  would  in  turn  cause  problems  during  the  transfer  to  G 

In  addition,  we  have  experimented  with  using  a  mixtnre  of  the  adaptive 
(r),  6)  strategy  with  the  non-adaptive  (p.q)  strategy  (cf.  section  3.2).  We 
have  found  a  (i|,q)  strategy  that  is  as  good  as  any  other  we  have  tried.  In 

this  strategy,  we  use  q  to  control  when  we  terminate  relaxation  on  a  certain 
grid  and  go  on  to  a  coarser  grid,  and  use  q  to  control  how  many  sweeps  to  do 
on  a  grid  after  transfer  from  a  coarser  grid  before  interpolating  onto  a  finer 
grid.  A  typical  set  of  parameters  that  worked  well  is  (q  =  .6,  q  =■  2).  The 
resulting  algorithm  is  fairly  insensitive  to  actual  values  of  q  and  q  and  is 
quite  robust.  It  is  also  quite  efficient.  It  consistently  achieved  an 
efficiency  of  less  than  about  12  units  per  order  of  magnitude  reduction  in  the 
residual  for  most  problems  that  we  have  encountered.  Some  of  these  problems 
have  very  singular  grids  which  presented  difficulties  for  all  of  the  previous 


algorithms. 
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7.  Summary 

In  this  paper,  we  study  arc-length  continnation  techniques  and  multi-grid 
techniques  for  solving  nonlinear  elliptic  eigenvalue  problems.  We  have 
applied  these  techniques  to  solve  a  model  nonlinear  elliptic  eigenvalue 
problem  (the  Bratu  problem).  We  have  found  that  as  long  as  we  stay  away  from 
singular  points,  the  two  techniques  combined  to  give  a  very  powerful  and 
efficient  procedure  for  tracing  solution  branches.  Near  singular  points, 
however,  the  standard  multi-grid  method  has  difficulty  converging  on  the 
linearized  elliptic  systems  that  arise  in  the  continuation  procedure.  One 
consequence  is  that  we  cannot  continue  past  the  limit  point  in  the  model 
problem.  This  divergence  is  successfully  analysed  and  several  modified 
multi-grid  algorithms  have  been  designed  based  on  this  analysis.  The  best  of 
these  modified  algorithms  performs  efficiently  and  reliably  arbitrarily  close 
to  the  singular  points.  This  enables  the  continuation  procedure  to  continue 
past  the  limit  point  with  no  difficulty.  It  seems  reasonable  that  this 
modified  multi-grid  algorithm  can  be  useful  in  more  general  situations  where 
nearly  singular  elliptic  systems  arise,  such  as  in  inverse  iteration  [11,  17]. 
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